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The statistical physics of homogeneous DNA is investigated by the imaginary time path integral 
formalism. The base pair stretchings are described by an ensemble of paths selected through a 
macroscopic constraint, the fulfillement of the second law of thermodynamics. The number of paths 
contributing to the partition function strongly increases around and above a specific temperature 
T* whereas the fraction of unbound base pairs grows continuosly around and above T* . The latter 
is identified with the denaturation temperature. Thus, the separation of the two complementary 
strands appears as a highly cooperative phenomenon displaying a smooth crossover versus T. The 
thermodynamical properties have been computed in a large temperature range by varying the size 
of the path ensemble at the lower bound of the range. No significant physical dependence on the 
system size has been envisaged. The entropy grows continuosly versus T while the specific heat 
displays a remarkable peak at T* . The location of the peak versus T varies with the stiffness of the 
anharmonic stacking interaction along the strand. The presented results suggest that denaturation 
in homogeneous DNA has the features of a second order phase transition. The method accounts for 
the cooperative behavior of a very large number of degrees of freedom while the computation time 
is kept within a reasonable limit. 



PACS numbers: 87.14.gk, 87.15.A-, 87.15.Zg, 05.10.-a 



I. Introduction 



Path integral methods provide a powerful tool to in- 
vestigate the properties of nonlinear dynamical systems 
with retarded interactions. Hamiltonian models with non 
local electron-phonon couplings, used in polymer science 
[l|, have been studied by such methods. In this paper the 
path integral approach [2| is applied to analyse the tem- 
perature driven DNA denaturation which occurs when 
the double stranded molecule separates into two coils 
while, at lower T, only localized openings exist |3|. Mod- 
elling of DNA melting has been motivated since long by 
the need to understand the transcription mechanism in 
which the double helix opens locally to allow a reading 
of the genetic code. In this regard, the Poland-Scheraga 
model |4| has been path-breakingand still makes the base 
of several investigations [1, H, 0, H, Ipl-IToll. Despite exten- 
sive analytical and numerical work jlll.ll2l. Il3l . [lil . [TEl [l6| , 
establishing character and nature of the melting phase 
transition, whether first or second order, remains a chal- 
lenging task [l7[. A significant advance towards a com- 
prehension of the DNA dynamics was made after Peyrard 
and Bishop (PB) introduced a ID Hamiltonian model [l8[ 
which recognized the role of nonlinearities in the molecule 
and reduced its great complexity to two essential 
interactions: i) a nonlinear coupling between the two 
bases in a pair connected by hydrogen bonds, ii) a har- 
monic stacking potential between adjacent bases along 
the strand. While solutions based on the transfer inte- 
gral method [201 ] predicted a smooth denaturation transi- 
tion both for homogeneous and heterogeneous DNA [2l[ , 
extensions of the PB model [13, [H| which include anhar- 
monic stacking interactions showed that the denaturation 
can be sharp, occurs at lower T than in the harmonic case 
and it is indeed a thermodynamic transition [24| . The 



latter statement does not contraddict general theorems 
on the impossibility of phase transitions in ID systems 
[25l ] as the PB Hamiltonian contains an on site potential 
depending on unbounded transverse stretchings. More- 
over, in the thermodynamic limit, the energy of a domain 
wall between two states of the molecule is infinite and an 
ordered state is then possible [26| . 

As the anharmonic PB model is assumed as fundamen- 
tal tool, it is understood that the emphasis of the present 
investigation is purely on modelling the unbinding of the 
two strands versus temperature. This is done by captur- 
ing the essentials of the interactions consistently with the 
spirit of the theory of critical phenomena. 

The PB Hamiltonian model is outlined in Section II 
focussing on those properties which make reliable to at- 
tack the problem by path integral techniques. Section 
III describes the method while the results concerning the 
character of the denaturation transition are discussed in 
Section IV. Some conclusions are drawn in Section V. 



II. Hamiltonian Model 

To begin, I consider a fragment of homogeneous DNA 
whose strands are represented by a set of point masses 
corresponding to the nucleotides. The Hamiltonian [22j 
for a chain of N base pairs (bps) with mass fi, reads: 
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V S (y n , Vn-l) = yff(j/n,J/n-l)(l/n ~ 2M-l) 2 
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where y n , the transverse stretching for the n- bp, mea- 
sures the relative pair separation from the ground state 
position and continuously grows from (closed bp) to 
oo (open bp after denaturation occurs). As the hydro- 
gen bond can be compressed with respect to the equilib- 
rium, y n takes also negative values with a lower cutoff 
accounting for the repulsive interaction of the phosphate 
groups. The longitudinal displacements are neglected in 
the model as their amplitudes are much smaller than the 
transverse stretchings 0, D and a, setting depth 

and width of the Morse potential VM(y n ), are site inde- 
pendent for homopolymer DNA. In the stacking poten- 
tial Vs(yn,yn-i), the factor g depends on the sum of 
the stretchings of two adjacent bases. K = pv 2 with 
v being the harmonic phonon frequency. For p = in 
Eq. (JTJ, the statistical mechanics can be worked out ex- 
actly through the transfer operator method [2(| which, 
in continuum approximation and strong K limit, defines 
a pseudo-Schrodinger equation for a particle in a Morse 
potential [28[. The latter has a discrete spectrum with 
localized eigenfunctions at T < T c being 
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the temperature at which the lowest eigenvalue van- 
ishes and even the ground state merges into the contin- 
uum. Kb is the Boltzmann constant. For T > T c only 
delocalized states exist. Thus, even the simplest har- 
monic stacking model shows that a transition occurs and 
T c is identified as the temperature at which the strands 
separate. Taking typical values for DNA parametriza- 
tion d, Hi] such as D = 30meV, a = 4.2l~ 1 and 

° —2 

K = 60meVA , one gets T c = 331K. Anharmonicity 
(p > 0) has a special meaning in DNA stacking. While 
the analytical form of g may not be unique in order to 
reproduce the denaturation transition [3(| , the key prop- 
erty of the stacking potential is the following: whenever 
either one of the bps is stretched over a distance larger 
than a -1 , the hydrogen bond breaks and the electronic 
distribution around the two pair mates is modified. Ac- 
cordingly, the stacking coupling (along each strand) be- 
tween neighboring bases in Eq. ([1]) drops from K(l+p) to 
K. Then, also the next bp tends to open as both bases are 
less closely packed along their respective strands. Here is 
the cooperative character which underlines the formation 
of a region with open bps, a bubble whose size increases 
with T. In heterogeneous DNA, also the length of the se- 
quence may affect the probability to open a bubble which 
is quantified in polymer network theories by the cooper- 
ativity parameter a: for sequences of intermediate length 
(up to ~ 10 4 bps), a small a (~ 10~ 4 — 10 -5 ) generally 
suppresses bubble formation while large portions of the 
helix unbind close to the melting making the transition 
highly cooperative [3l|, HH, [33|. For long sequences, a 
is larger and small bubbles may form already below the 
transition which accordingly is expected to be less coop- 
erative. 



When the stacking decreases the vibrational mode be- 
tween the two bases softens thus reducing its contribu- 
tion to the free energy. Then, as confirmed by molecu- 
lar dynamics simulations [28| . the denaturation onset is 
signalled by a phonon mode softening which, in general 
[34 , [35|], may point to the occurence of a phase tran- 
sition. For these reasons, anharmonic stacking models 
well account for cooperativity effects in the formation of 
open domains. Being the latter a collective phenomenon, 
computational methods have to include hundreds of base 
pairs in order to study the dynamics even of a molecule 
fragment and this requires considerable computational 
power. These facts suggest that models at intermediate 
scales are important and motivate the idea to investigate 
the DNA denaturation by path integral techniques. 



III. Path Integral Method 

The path integral formalism defines the time (t) evolu- 
tion amplitude between two points, say "a" and "b", as 
a sum over all histories along which a system can evolve 
in going from "a" to "b". Each history is weighed by 
a phase factor, the exponential of the action associated 
to a given path 36]. At finite temperature, the thermal 
properties of the system can be derived by weighing the 
contributions by the particle paths x(t) running along 
the imaginary time axis r, after an analytic continua- 
tion is performed: r = it. Accordingly, in the statistical 
formalism, the Euclidean action A[x(t)] replaces the me- 
chanical canonical action and the partition function is an 
integral in the path phase space. Each path is weighed 
by a probability factor exp(— A[x(r)]). Only closed paths 
contribute to the statistical partition function, the inte- 
gration being a trace integration [53] . In the calculations, 
not all histories can be accounted for and, given the spe- 
cific problem, one has to select the suitable class of paths 
which mainly contribute (are expected to contribute) to 
the physical properties. 

For the model in Eq. ([T]) , the path integral description 
naturally follows by replacing the bps transverse stretch- 
ing by a one dimensional path x which is continuous func- 
tion of r. Accordingly, the space interaction are mapped 
onto the time scale as 
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x(t), y„_x — ► x(t') , t' = t - At , (3) 



with r G [0, 0] and f3 being the inverse temperature. 
While only adjacent bases are taken in Eq. ([1]), the stack- 
ing interactions range can be varied in the path integral 
model by tuning the retardation At in Eq. (J3]) . The dis- 
crete lattice nature of the Hamiltonian is maintained by 
the replacement in Eq. ([3]). In fact, for any (3, the com- 
putation is made convergent by taking a large but finite 
number N T of points in the r range. The path x(t) is 
further expanded in Fourier series with cutoff Mp 
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x(t) = Xq + 



Mf 

E 



m cos(w m r) + b m sin(w m T) (4) 



and u> m — 2rmr/f3. As x(0) = x(/3), periodic bound- 
ary conditions hold for the present method analogously 
to those imposed for the ID finite chain described by 
Eq. ([T]). The heart of the matter lies in summing over a 
suitable set of paths such that the computation of the 
thermodynamical properties becomes sufficiently accu- 
rate. As explained hereafter, the selection of the paths 
ensemble mainly contributing to the partition function 
is done on the base of physical constraints. By Eq. ([3]), 
the index n maps onto a single r value and, by Eq. ([J), 
for any r an ensemble of path coefficients is taken in the 
computation. 

The imaginary time partition function [2| is given by 



Z = I Sx(t) exp[-A(x(r))] 



A(x(t)) 



dr 



^x{rf+V{x{r)) 
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»*(r) = ^- / dx n (^) 2 / da m J db m , 

(5) 



where A(x(t)) is the Euclidean action for a particle in 
the potential V(x(t)). Dx(t) is the measure of integra- 
tion which normalizes the free particle action 



Dx(t) exp 



dr— ±(r) 2 
2 v ' 



= 1 



(6) 



and § denotes integration over closed particle trajecto- 
ries consistently with periodic boundary conditions hold- 
ing for paths in Eq. (j4]). 

A M is the thermal wavelength whose form depends 
on the model whether quantum or classical. DNA de- 
naturation occurs at temperatures for which a classical 
model applies. Then, the time derivative y n (Eq. ([1])) 
maps onto the imaginary time derivative x(r) (Eq. l[5])). 
the proper replacement being: d/dt — > {v(3)d/dr hence, 
Aju = y/n/(3K. Note that the pseudo-Schrodinger equa- 
tion is also solved in a classical framework with h replaced 
by { V f3)-\ 

Applying Eq. ([3]) to the Hamiltonian in Eq. ((T|), Z in 
Eq. ^ transforms into 



Z = 



Uo M F 

dx Yl 



ds r 



dt r . 



exp(-s^ - t 2 n )E(x , s m ,t m ) 



E(x ,s m ,t m ) = exp^-y drV^ T {x) 

Vat(x) = Vm{x) + V s {x,x') 
Vm(x) = D[exp(-aar) - l] 2 

—g[x,x\{x 



V s (x,x') 

g[x,x']= l + pexp[— a 
x = x(t) ; x' = x(t') 
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where the measure has been written in a form suitable 
for computation by introducing the cutoffs Uq and U on 
the Fourier coefficients integrals [39]. By inserting the 
expansion in Eq. ([!]), the l.h.s. of Eq. ^ transforms 
into a product of Gaufiian integrals. Normalization of 
the latter gives the mathematical criterion to set Uq and 
U through the conditions: 



U 



dxo = 1 



-U 



ds m exp(-s; 



(8) 



While Uq is set by the first in Eq. ([8]), U can be deter- 
mined after a series expansion for the Gaufiian integral 
[HI is applied to the l.h.s. in the second of Eq. ((8]). How- 
ever, there is no a priori reason why the mathematical 
cutoffs had to fulfill also the physical requirements of the 
problem in Eq. (JT]). In fact, as extensively discussed be- 
low, too negative Fourier coefficients would induce large 
negative path amplitudes and diverging Vm which, in 
turn, would yield zero contribution to Z. Such paths, 
although allowed by Eq. ([S]), have therefore to be dis- 
carded as their inclusion in the computation would pro- 
duce a decreasing entropy versus T. While this event 
cannot be accepted on general physical grounds, it also 
provides clue to solving consistently the cutoff problem: 
imposing the second law of thermodynamics, I solve the 
path integral in Eq. ((7|) selecting at any T the suitable 
path ensemble for the transverse stretchings. Let's see in 
detail the technicalities of the method. 

Besides the five model potential parameters, namely 
D, a, K, p, a, the path integral method contains some 
intrinsic input parameters. The latter are: a) The cutoff 
Mp. b) The retardation At. c) The number of points 
N T for the dr-integral in Eq. 0. d) The number of 
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integration points over the Fourier coefficients xq, s m and 

tra- 
gi.) Mp = 1 suffices to make the calculation conver- 
gent in the sense that Z would not change by taking 
Mf = 2. This result can be understood as follows: the 
inclusion of higher Fourier components leads to larger 
(absolute) values for the path amplitudes but Vm takes 
care of suppressing paths that leave the reasonable bound- 
aries. While there is no reason to include paths whose 
fate is that of being suppressed, a proper ensemble of 
paths can be built by taking a sufficiently high number 
of integration points for the first Fourier component in 
a sufficiently broad range. This ensures that the pro- 
gramme samples a large portion of the phase space in 
which all the relevant path amplitudes are accounted for. 
Incidentally, the (Mp — l)-strategy also permits to save 
a great amount of computing time. 

b) This investigation is restricted to adjacent bps inter- 
actions along the strand. Accordingly, r and t' in Eq. ([3]), 
are first neighbors in the discrete imaginary time lattice 
and the retardation is: At = (3/N T . 

c) N T = 100 suffices to make the action numerically 
stable. This value also sets the number of bps as N T co- 
incides with N in Eq. (TTJ) due to the mapping in Eq. ([3]). 
While N — 100 is a (minimal) reliable choice for a co- 
operative phenomenon such as denaturation, it turns out 
that larger N T (N) would not change anything in the 
computation of Eq. |(7J). At this stage, the reader may 
wonder whether and how the path integral method can 
account for those finite size effects [53, |41j which are con- 
sidered relevant in studies of the DNA properties. The 
fact is that, for the present method, N is not the exclu- 
sive tuning parameter in order to set the system size, the 
latter being rather determined by the total number of 
paths allowed for the N-bps. Thus, although N is kept 
constant, the system size may still increase/decrease by 
considering a larger/smaller path ensemble which even- 
tually works as the true scaling parameter in the path 
integral method. Here, we see the deep meaning of the 
method under lied by the mapping in Eq. ((3]). The n — bp 
maps onto an effective ensemble of paths N e f f taken at a 
specific r. The latter may be viewed in itself as an ensem- 
ble of base pairs. As this holds for any of the N T points, 
the true scaling parameter turns out to be N T ■ N e ff. 

d) Taking one Fourier component, the program eval- 
uates at any T the contribution to Z by three integrals 
over Fourier coefficients in Eq. ([7])- If, for a given number 
of integration points, the increasing entropy constraint is 
fulfilled then the good paths are included in the compu- 
tation. The number of good paths for a single r is N e ff. 
The total number of good paths over which computations 
are carried out is N T ■ N e f / . 

Clearly, the computing time depends on the num- 
ber of paths required to stabilize the system. Many 
paths should contribute to the thermodynamical proper- 
ties consistently with the previous observations regard- 
ing the cooperativity behavior in the double helix. Such 
number is expected to grow markedly around and above 
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FIG. 1: (Color online) (a) Lower and upper cutoffs (in A) for 
the paths, representing the base pairs transverse stretchings, 
versus Temperature ; (b) Morse potential Vm (in meV) versus 
path amplitude x. D is in meV and a in A . 



the transition. This would be a consistent hallmark for 
the reliability of the path integral approach to the DNA 
Hamiltonian. 



IV. Discussion 

The issue of the maximum (absolute) path amplitudes 
is key to understand the double helix denaturation [42j 
and it has to be discussed consistently with the available 
experimental data. As put forward above, the properties 
of the potential Vm, plotted in Fig. QJb), suggest which 
path ensemble is physically meaningful. 

Hydrogen bond stretchings are generally broken for 
lengths ~ 2A. For such separation, the potential is flat 
with a plateau value D ~ 30meV and the force between 
the bases vanishes. Accordingly, larger amplitudes would 
yield no contribution. Then, the bp breaking energy D is 
~ KbT above room temperature. This justifies the pa- 
rameter values used in Fig. [TJb) . Largely positive paths 
would sample the flat portion of Vm and add constant 
terms to the action which do not affect the free energy 
derivatives. 

On the other side, steric hindrance processes prevent 
x(t) from being too negative. Paths x(t) ~ — 0.14A 
also experience a value D ~ SOmeV while larger negative 
paths encounter the hard core repulsive side of Vm ■ Here 
the latter grows exponentially thus yielding a large action 
which, in turn, makes a vanishing contribution to Z . 

These are the intuitive reasons to set a finite range 
which inspire the computational method described in 
the previous Section. Moreover, it is physically plau- 
sible that such range had to be modulated by tem- 
perature effects allowing for larger path amplitudes at 
higher T. Saying T* the denaturation temperature, 
the path displacements may reasonably vary in the 
range x(r) e [x min , x max ] with x min ~ -0.14A and 
Xmax ~ 2A around the crossover. The functions: 
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FIG. 2: (Color online) (a) Effective Number of Paths con- 
tributing to the partition function for three stacking parame- 
ters, a in A . (b) Specific heat computed via Eq. 0. 



Xmin = -0.14exp[(T-T*)/T*] and x max = 2exp[(T- 
T*)/T*\ , provide suitable (although not unique) expres- 
sions to account for linear temperature dependence in a 
window around T* whereas larger deviations occur out- 
side that window. x m i n and x max are plotted in Fig.fija) 
taking for T* the value given by Eq. ©. However, this 
assumption does not imply that the melting occurs pre- 
cisely at T c . In fact the true denaturation temperature 
(found below) does not coincide with the above given T c 
which stems from an analytic harmonic model. More- 
over, the specific choice for upper and lower bounds has 
not to imply in principle that a transition exists nor it 
has to force the transition to take place at T* . Certainly, 
if the denaturation occurs, a significant fraction of bps 
is expected to be ~ 2A. This feature will be hereafter 
investigated. 

Thus, after setting a temperature value, the program 
searches for those paths {x(t),x(t')} which lie inside the 
range [x m i n , x max \ for any r € [0,/3] and adds their con- 
tribution to Z . Such a check is done for any set of Fourier 
coefficients which in fact defines a specific path. Had the 
entropy not to have a positive derivative, the computa- 
tion should be re-run for a new number of integration 
points which, in turn, leads to re-build the ensemble of 
Fourier coefficients and to repeat the procedure to fil- 
ter the suitable paths. It is seen numerically that the 



restriction due to x max can be lifted without affecting 
the physical properties whereas the lower bound x m i n is 
required in the computational method. 

As a final output N e f f is obtained versus T as shown 
in Fig. Uta) for the same values D, a, K as in Figs. [TJ 

and for three stacking parameters p with a = 0.35A 

While all plots begin with N eff ~ 8000 at T = 
250K, their slopes display a remarkable dependence on 
p and, as a general feature, N e ff is higher for stronger 
anharmonicities. In particular, for the case p = 2, N e ff 
increases at T ~ 315 A. This feature is mirrored in the 
plot (Fig.[^b)) of the specific heat which is a significant 
hallmark for denaturation as it is proportional to the 
differential melting curve (43|. For p — 2, a peak is 
found at the same T and in general, the specific heat 
displays a sharp increase when the degree of cooperativity 
(measured by N e ft) markedly grows. The peak shifts 
towards lower T by increasing p. This finding, shared by 
previous studies [22j , is physically meaningful as systems 
with larger anharmonic couplings should have a higher 
number of pairs which break at lower T. 



A. Phase Transition 

While the results displayed so far witness by them- 
selves reliability and consistency of the method, they also 
induce to face the unresolved issue regarding the order 
of the denaturation transition. Two questions need to be 
addressed: 

1) how does the entropy behave? 

2) To which extent do the thermodynamical properties 
depend on the number of paths included in the calcula- 
tion? Or, how does the macroscopic system scales versus 

Taking the intermediate plot, p = 1 in Fig. [2l I have 
thus progressively increased N e f t over the value N e f t = 
9920 at T = 260K in Fig. [^a). While the T axis in the 
latter starts from 250K, the 1QK shift upwards allows one 
to distinguish among the three plots. T = 260-ftT is then 
the lower bound in the T window considered hereafter. 
Once the starting N e ff is set at the lower bound, the 
thermodynamical properties are evaluated at higher T 
according to the procedure detailed above. 

In Fig. [31 the calculation begins with N e ff — 31792, 
a value larger by a factor three than in Fig. O Fig. [H(a) 
plots: i) the total number of paths versus T starting with 
~ 3.2 • 10 6 at T — 260-R"; ii) the number of paths whose 
amplitude is larger than 1A; Hi) the number of paths 
whose amplitude is larger than 2A. All curves increase 
versus T and show a kink at T*. This feature is more 
pronounced and evident than in Fig. [2] confirming that 
denaturation is indeed a highly cooperative phenomenon. 
While paths > 1A already sample the plateau of Vm 
where the bps unbind, paths > 2A certainly belong to 
broken pairs. Accordingly, the latter become more nu- 
merous around and above T*. Thus, approaching the 
transition, the system incorporates a larger number of 
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paths with an increasing fraction of broad path ampli- 
tudes. As such fractions are considered reliable indicators 
of the order of the denaturation transition, quantitative 
estimates are greatly relevant. As the inset makes evi- 
dent, the fractions of paths > 1A and > 2A respectively 
(normalized over N T ■ N e ff), grow continuosly versus T. 
Around T*, one third of all paths are larger than lA: it 
means that either they are broken or about to break. We 
see no step-like increase at T* and this behavior is fully 
consistent with the continuous entropy growth shown in 
Fig. [3jb) . The slight irregularity appearing in the en- 
tropy plot at T* is responsible for the remarkable peak 
in the specific heat plotted in Fig. [J2c) . The peak remains 
pinned at T* = 363K as in Fig. H(b). The inset shows 
an enlargement with a temperature resolution which ap- 
preciates O.IK. Finer partitions may be taken further 
increasing computing time. 

The results presented so far point to the occurence of a 
smooth crossover at T*. While evidence is emerging for 
classyfing the transition as of second order, the careful 
reader may argue that the aforementioned question 2) 
is not yet fully answered. I have then further increased 
the starting N e f f : all the obtained curves look like those 
presented in Fig. [3] with, at most, shifts of a few degrees 
in the location of T* . Fig. [4] shows the results obtained 
for N e ff = 47494 at T = 260K: in this case T* = 352K 
which makes the largest observed shift with respect to 
Figs. [Hand [3] Note that ~ 8 • 10 6 paths participate to 
the computation around T*, with 37% being larger than 
lA. Remarkably, ~ 18 • 10 6 paths are accounted for at 
the largest T here considered. 

Some features shown in the plots and drawn from the 
whole numerical work are in order: 

*) The entropy is always a continuosly increasing func- 
tion of T whatever the starting value for N e ff may be. 

*) At T*, the entropy always displays a slightly more 
pronounced enhancement although there is no finite melt- 
ing entropy which would point to a first order phase tran- 
sition [2J] . The small entropy kink is instead responsible 
for the peak in the specific heat. The crossover is smooth, 
the phase transition appears to be of second order. 

*) The entropy values are larger by including more 
paths in the computation. 

*) Above T*, the entropy goes on growing but with 
smaller gradient consistently with the fact that, as the 
strands separation proceeds, the probability to add en- 
tropic effects is reduced. 

*) The location of T* slightly varies with increasing 
N e ff but the changes are not significant oscillating in 
the range T* ~ [350 - 365K] for p = 1. 

While different choices for the stiffness parameters in 
the PB Hamiltonian are not expected to change these 
qualitative physical features, it may be instead worth in- 
vestigating by path integrals also Hamiltonian models 
with different stacking potentials both for homogeneous 
and heterogeneous DNA [441 ]. 
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FIG. 3: (Color online) Anharmonic stacking case p = 1. 
N e ff — 31792 at T = 260K is assumed. T axis starts with a 
10K shift with respect to Fig. [2]in order to distinguish among 
the three plots corresponding to different p values, (a) Paths 
larger than 1A, 2A and Total Number of paths (N T -N e ff) are 
plotted. Inset: ratios of paths whose amplitude is larger than 
lA and 2A respectively, over N T ■ N e ff. (b) Entropy versus 
Temperature, (c) Specific Heat versus Temperature. Inset: 
Magnification of the Specific Heat Peak with a temperature 
resolution of 0.1K. 



V. Conclusion 

The homogeneous DNA denaturation has been studied 
in the path integral formalism by mapping the essential 
interactions of the Peyrard-Bishop Hamiltonian onto the 
imaginary time scale. While the nonlinearities in the 
stacking interactions are recognized as a key feature of 
the Hamiltonian, the path integral method permits to 
deal with the highly cooperative character of the denat- 
uration by incorporating a very high number of degrees 
of freedom. I have thought of the transverse stretchings 
for the base pairs as time dependent paths: in the imag- 
inary time formalism, this amounts to incorporate the 
temperature effects in the base pairs displacements. The 
path amplitudes are taken consistently with the shape 
of the Morse potential modelling the hydrogen bonds. 
Accordingly, I have set up a computational method to 
monitor, versus temperature, the paths ensemble which 
mainly contributes to the partition function, hence to 
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FIG. 4: (Color online) As in Fig. [3] but with N eff = 47494 
at T = 260 AT. 

the macroscopic equilibrium properties of the system. A 
temperature window has been considered in which de- 
naturation is expected to occur. Given the size of the 
paths ensemble at the lower bound of such window, the 
second law of thermodynamics has been invoked to build 
the ensemble at any larger T. The method works self- 
consistently selecting the minimum number of paths re- 
quired (at any T) to pursue the goal, that is the growing 
entropy constraint. There is no need to add further paths 



to the partition function once such goal is achieved: I feel 
that a minimum effort principle should hold. Instead, 
the system physical properties may depend on the initial 
ensemble size and care has been taken to unravel this is- 
sue. Thus, the partition function and its derivatives have 
been computed by varying the boundary condition but 
no significant variation in the physical behavior has been 
observed. It is understood that the initial size should 
be large enough to allow for the appearance of coop- 
erative effects. Remarkably, the ensemble size strongly 
grows around and above the denaturation temperature. 
The denaturation of the complementary strands always 
shows up as a smooth crossover whose location along 
the T-axis may depend on the stiffness of the stacking 
interaction. Focussing on the system with moderate an- 
harmonic stacking, it has been found that: i) the entropy 
grows continuosly versus T , ii) the specific heat displays 
a peak precisely at the same temperature for which the 
path ensemble size sharply grows. 

These results suggest that the denaturation in homo- 
geneous DNA is a highly cooperative phenomenon with 
the hallmarks of a second order phase transition. Consis- 
tently, the fraction of paths sampling the plateau of the 
Morse potential grows continuosly also around the de- 
naturation temperature. No step-like increase has been 
observed in the paths fractions which describe the un- 
binding of the base pairs. Eventually, I wish to empha- 
size that the thermodynamical constraint applied on the 
system does not force a priori the entropy to grow con- 
tinuosly. While the system autonomously selects at any 
T the total number of paths to fulfill such constraint, 
both a step discontinuity and a continuous behavior may 
in principle appear. The obtained plots are therefore not 
an artifact of the method. 

After testing the feasibility of the path integral ap- 
proach to a large size system as homogeneous DNA, I 
feel that the presented formalism may be further im- 
proved/extended to account for bubble formations in het- 
erogeneous models with inhomogeneous sequences. 
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